cg_regions <- read_rds(here::here("data/derived_data/cg_regions.rds"))
animal_regions <- read_rds(here::here("data/derived_data/animal_regions.rds"))
#prism_zip <- read_rds(here::here("data/derived_data/prism_zip.rds"))
cg_regions %>%
filter(!region %in% c(4, 6)) %>%
filter(trait == "ww" & var == "cg_sol") %>%
filter(n_animals > 4) %>%
group_by(year) %>%
mutate(year_sd = sd(value),
year_mean = mean(value),
keep = case_when(
value > year_mean + year_sd ~ "YES",
value < year_mean - year_sd ~ "YES",
TRUE ~ "NO"
)) %>%
ungroup() %>%
filter(keep == "YES") %>%
# sample_frac(0.2) %>%
ggplot(aes(
x = lng,
y = lat,
color = "goldenrod",
size = n_animals
)) +
geom_point(
alpha = 0.3,
#size = 0.5
) +
scale_size(
range = c(1,4)
) +
usa +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
coord_map("albers", lat0 = 39, lat1 = 45) +
cowplot::theme_map() +
labs(
x = NULL,
y = NULL,
color = NULL,
title = str_wrap( "Weaning weight CG solutions, low to high", width = 55)
) +
#Set the "anchoring point" of the legend (bottom-left is 0,0; top-right is 1,1)
#Put bottom-left corner of legend box in bottom-left corner of graph
theme(
axis.title.x = element_blank(),
axis.title.y = element_blank(),
# Padding around the legend
# top, right, bottom, left
# legend.box.margin = margin(b = 0, r = 0.2, unit = "cm"),
# Padding around the plot
plot.margin = margin(
t = 0.7,
r = 0,
b = 0.7,
l = 1,
unit = "cm"
),
plot.title = element_text(
# size = 56,
# family = "lato",
size = 22,
vjust = 6,
face = "italic"
)
) +
# https://stackoverflow.com/questions/32656553/plot-legend-below-the-graphs-and-legend-title-above-the-legend-in-ggplot2
guides(
color = FALSE,
size = FALSE
) +
gganimate::transition_reveal(value)
##
Frame 1 (1%)
Frame 2 (2%)
Frame 3 (3%)
Frame 4 (4%)
Frame 5 (5%)
Frame 6 (6%)
Frame 7 (7%)
Frame 8 (8%)
Frame 9 (9%)
Frame 10 (10%)
Frame 11 (11%)
Frame 12 (12%)
Frame 13 (13%)
Frame 14 (14%)
Frame 15 (15%)
Frame 16 (16%)
Frame 17 (17%)
Frame 18 (18%)
Frame 19 (19%)
Frame 20 (20%)
Frame 21 (21%)
Frame 22 (22%)
Frame 23 (23%)
Frame 24 (24%)
Frame 25 (25%)
Frame 26 (26%)
Frame 27 (27%)
Frame 28 (28%)
Frame 29 (29%)
Frame 30 (30%)
Frame 31 (31%)
Frame 32 (32%)
Frame 33 (33%)
Frame 34 (34%)
Frame 35 (35%)
Frame 36 (36%)
Frame 37 (37%)
Frame 38 (38%)
Frame 39 (39%)
Frame 40 (40%)
Frame 41 (41%)
Frame 42 (42%)
Frame 43 (43%)
Frame 44 (44%)
Frame 45 (45%)
Frame 46 (46%)
Frame 47 (47%)
Frame 48 (48%)
Frame 49 (49%)
Frame 50 (50%)
Frame 51 (51%)
Frame 52 (52%)
Frame 53 (53%)
Frame 54 (54%)
Frame 55 (55%)
Frame 56 (56%)
Frame 57 (57%)
Frame 58 (58%)
Frame 59 (59%)
Frame 60 (60%)
Frame 61 (61%)
Frame 62 (62%)
Frame 63 (63%)
Frame 64 (64%)
Frame 65 (65%)
Frame 66 (66%)
Frame 67 (67%)
Frame 68 (68%)
Frame 69 (69%)
Frame 70 (70%)
Frame 71 (71%)
Frame 72 (72%)
Frame 73 (73%)
Frame 74 (74%)
Frame 75 (75%)
Frame 76 (76%)
Frame 77 (77%)
Frame 78 (78%)
Frame 79 (79%)
Frame 80 (80%)
Frame 81 (81%)
Frame 82 (82%)
Frame 83 (83%)
Frame 84 (84%)
Frame 85 (85%)
Frame 86 (86%)
Frame 87 (87%)
Frame 88 (88%)
Frame 89 (89%)
Frame 90 (90%)
Frame 91 (91%)
Frame 92 (92%)
Frame 93 (93%)
Frame 94 (94%)
Frame 95 (95%)
Frame 96 (96%)
Frame 97 (97%)
Frame 98 (98%)
Frame 99 (99%)
Frame 100 (100%)
## Finalizing encoding... done!
# anim_save(filename = here::here("figures/ww_solutions_appear.gif"), nframes = 48, fps = 2)
cg_regions %>%
filter(!region %in% c(4, 6)) %>%
solutions_line(
effect_var = "cg_sol",
trait_var = "ww",
y_lab = "Mean CG solution",
plot_title = "Weaning weight contemporary group solutions reflect year-to-year environmental trends"
) +
lims(y = c(440, 580))
#ggsave(here::here("figures/ww_cg.solutions_line_year.png"), width = 10, height = 5.4)
animal_regions %>%
filter(!region %in% c(4, 6)) %>%
solutions_line(
effect_var = "bv_sol",
trait_var = "ww",
y_lab = "Mean breeding value",
plot_title = "Weaning weight genetic trend: 1973-2019"
)
#ggsave(here::here("figures/ww_bv.solutions_line_year.png"), width = 10, height = 5.4)
animal_regions %>%
filter(!region %in% c(4, 6)) %>%
solutions_line(
effect_var = "mat_sol",
trait_var = "ww",
y_lab = "Mean maternal effect solution",
plot_title = "Milk genetic trend: 1973-2019"
)
#ggsave(here::here("figures/ww_mat.solutions_line_year.png"), width = 10, height = 5.4)
(Removed data from 1972: weird outlier, only one region with data)
cg_regions %>%
filter(!region %in% c(4, 6)) %>%
solutions_line(
effect_var = "cg_sol",
trait_var = "pwg",
y_lab = "Mean CG solution",
plot_title = "Post-weaning gain contemporary group solutions reflect year-to-year environmental trends"
)
#ggsave(here::here("figures/pwg_cg.solutions_line_year.png"), width = 10, height = 5.4)
animal_regions %>%
filter(!region %in% c(4, 6)) %>%
solutions_line(
effect_var = "bv_sol",
trait_var = "pwg",
y_lab = "Mean breeing value",
plot_title = "Post-weaning gain genetic trend: 1973-2019"
)
#ggsave(here::here("figures/pwg_bv.solutions_line_year.png"), width = 10, height = 5.4)
cg_regions %>%
filter(!region %in% c(4, 6)) %>%
filter(trait == "pwg" & var == "cg_sol") %>%
mutate(yr = lubridate::year(weigh_date),
region = as.character(region)) %>%
group_by(region, yr) %>%
summarise(mean = mean(value)) %>%
bind_rows(
cg_regions %>%
filter(trait == "pwg" & var == "cg_sol") %>%
mutate(yr = lubridate::year(weigh_date)) %>%
group_by(yr) %>%
summarise(mean = mean(value)) %>%
mutate(region = "All regions")
) %>%
arrange(yr)
cg_regions %>%
filter(!region %in% c(4, 6)) %>%
solutions_var_line(
effect_var = "cg_sol",
trait_var = "ww",
plot_title = "The variance of weaning weight CG solutions changes little across regions or time"
)
#ggsave(here::here("figures/ww_cg.solutions_variance.png"), width = 10, height = 12)
cg_regions %>%
filter(!region %in% c(4, 6)) %>%
solutions_var_line(
effect_var = "cg_sol",
trait_var = "pwg",
plot_title = "The variance of post-weaning gain CG solutions changes little across regions or time"
)
#ggsave(here::here("figures/pwg_cg.solutions_variance.png"), width = 10, height = 12)
p1 <-
cg_regions %>%
filter(!region %in% c(4, 6)) %>%
filter(trait == "ww" & var == "cg_sol") %>%
filter(n_animals > 4) %>%
mutate(yr = lubridate::year(weigh_date),
yr = as.integer(yr)) %>%
group_by(region) %>%
mutate(comp = mean(value)) %>%
ungroup() %>%
group_by(zip, yr) %>%
mutate(sz = sum(n_animals),
zip_sol = mean(value)) %>%
ungroup() %>%
mutate(
overunder =
case_when(
zip_sol > comp ~ "Above mean",
zip_sol < comp ~ "Below mean"
),
howmuch = zip_sol - comp,
howmuch = abs(howmuch)
) %>%
select(trait, var, yr, region, comp, zip, sz, overunder, howmuch, lng, lat) %>%
distinct() %>%
filter(howmuch < 300) %>%
solutions_map(
effect_var = "cg_sol",
trait_var = "ww",
plot_title = "test",
color_var = howmuch,
size_var = sz,
size_range = c(1, 4)
) +
facet_wrap(~ overunder) +
theme(legend.position = "bottom",
strip.text.x = element_text(size = 16)) +
gganimate::transition_time(yr) +
labs(title = "Weaning weigh CG solutions: {frame_time}")
anim_save(filename = here::here("figures/ww_overunder_byregion.gif"), animation = plot, nframes = 48, fps = 1)
p2 <-
cg_regions %>%
filter(!region %in% c(4, 6)) %>%
filter(trait == "ww" & var == "cg_sol") %>%
mutate(yr = lubridate::year(weigh_date),
yr = as.integer(yr)) %>%
group_by(yr) %>%
mutate(comp = mean(solution)) %>%
ungroup() %>%
group_by(zip, yr) %>%
mutate(sz = sum(n_animals),
zip_sol = mean(solution)) %>%
ungroup() %>%
mutate(
overunder =
case_when(
zip_sol > comp ~ "Above mean",
zip_sol < comp ~ "Below mean"
),
howmuch = zip_sol - comp,
howmuch = abs(howmuch)
) %>%
select(trait, var, yr, region, comp, zip, sz, overunder, howmuch, lng, lat) %>%
distinct() %>%
filter(howmuch < 300) %>%
solutions_map(
effect_var = "cg_sol",
trait_var = "ww",
plot_title = "test",
color_var = howmuch,
size_var = sz,
size_range = c(1, 4)
) +
facet_wrap(~ overunder) +
theme(legend.position = "bottom",
strip.text.x = element_text(size = 16)) +
gganimate::transition_time(yr) +
labs(title = "Weaning weigh CG solutions: {frame_time}")
#anim_save(filename = here::here("figures/ww_overunder_byyear.gif"), animation = plot_byyear, nframes = 48, fps = 1)
cg_regions %>%
filter(var == "cg_sol" & trait == "ww") %>%
mutate(year = lubridate::year(weigh_date)) %>%
left_join(prism_zip %>%
select(-lat,-lng) %>%
mutate(year = year - 1)) %>%
group_by(year) %>%
mutate(
yr_sol = mean(value),
yr_temp = mean(tmean),
yr_rain = mean(ppt),
sd_sol = sd(value),
sd_temp = sd(tmean),
sd_rain = sd(ppt),
) %>%
ungroup() %>%
filter(region %in% c(1, 2, 5)) %>%
group_by(region, year) %>%
summarise(
dev_sol = (mean(value) - unique(yr_sol))/unique(sd_sol),
dev_temp = (mean(tmean) - unique(yr_temp))/unique(sd_temp),
dev_rain = (mean(ppt) - unique(yr_rain))/unique(sd_rain)
) %>%
ungroup() %>%
reshape2::melt(id = c("region", "year")) %>%
ggplot(
aes(
x = year,
y = value,
linetype = variable,
color = forcats::as_factor(region)
)
) +
geom_line(size = 1.25) +
scale_linetype_manual(
values = c(
"dev_sol" = "solid",
"dev_temp" = "dotted",
"dev_rain" = "twodash"
),
labels = c(
"dev_sol" = "CG solution",
"dev_temp" = "Temperature",
"dev_rain" = "Rain"
)
) +
scale_color_manual(
values = c(
"1" = "tomato2",
"2" = "darkslategray4",
"3" = "springgreen3",
"4" = "brown",
"5" = "goldenrod1",
"6" = "gray50",
"7" = "deeppink3",
"8" = "gray17",
"9" = "slateblue2"
),
labels = c(
"1" = "1: Desert",
"2" = "2: Southeast",
"3" = "3: High Plains",
"4" = "4: Rainforest",
"5" = "5: Arid Prairie",
"6" = "6: Cold Desert",
"7" = "7: Forested Mountains",
"8" = "8: Fescue Belt",
"9" = "9: Upper Midwest & Northeast"
)
) +
theme_classic() +
theme(
plot.title = element_text(
size = 22,
face = "italic",
margin = margin(
t = 0,
r = 0,
b = 13,
l = 0
)
),
axis.title = element_text(
size = 16),
axis.title.y = element_text(margin = margin(
t = 0,
r = 13,
b = 0,
l = 0
)),
axis.title.x = element_text(margin = margin(
t = 13,
r = 0,
b = 0,
l = 0
)),
axis.text = element_text(
size = 14),
legend.text = element_text(
size = 14)
) +
labs(x = NULL,
y = "Standardized value",
color = NULL,
linetype = NULL,
title = str_wrap("Environmental variables are not directly related to weaning weight CG solutions", width = 55)
) +
facet_wrap(~ region,
nrow = 3)
cg_regions %>%
filter(var == "cg_sol" & trait == "ww") %>%
mutate(year = lubridate::year(weigh_date),
) %>%
left_join(prism_zip %>%
select(-lat, -lng) %>%
mutate(year = year - 1)) %>%
group_by(region, desc, year) %>%
summarise(
yr_sol = mean(value),
yr_temp = mean(tmean),
yr_rain = mean(ppt),
sd_sol = sd(value),
sd_temp = sd(tmean),
sd_rain = sd(ppt),
) %>%
filter(region %in% c(1, 2, 5)) %>%
ungroup() %>%
mutate(
desc = str_remove(desc, "\\([[:alpha:]]+\\)")
) %>%
ggplot(aes(
y = yr_temp,
x = yr_sol,
color = forcats::as_factor(region))) +
geom_point() +
geom_smooth(method = "lm") +
scale_color_manual(
values = c(
"1" = "tomato2",
"2" = "darkslategray4",
"3" = "springgreen3",
"4" = "brown",
"5" = "goldenrod1",
"6" = "gray50",
"7" = "deeppink3",
"8" = "gray17",
"9" = "slateblue2"
),
labels = c(
"1" = "1: Desert",
"2" = "2: Southeast",
"3" = "3: High Plains",
"4" = "4: Rainforest",
"5" = "5: Arid Prairie",
"6" = "6: Cold Desert",
"7" = "7: Forested Mountains",
"8" = "8: Fescue Belt",
"9" = "9: Upper Midwest & Northeast"
)
) +
theme_classic() +
facet_wrap(~ desc,
nrow = 3,
scales = "free_y") +
theme(
plot.title = element_text(
size = 22,
face = "italic",
margin = margin(
t = 0,
r = 0,
b = 13,
l = 0
)
),
axis.title = element_text(
size = 16),
axis.title.y = element_text(margin = margin(
t = 0,
r = 13,
b = 0,
l = 0
)),
axis.title.x = element_text(margin = margin(
t = 13,
r = 0,
b = 0,
l = 0
)),
axis.text = element_text(
size = 16),
legend.text = element_text(
size = 14),
strip.text.x = element_text(size = 14)
) +
labs(x = "Mean CG solution for the year",
color = NULL,
y = "Mean temperature for the year",
title = str_wrap("Mean annual temperature is not directly related to mean annual weaning weight CG solution",
width = 38))